Contents

Land Cover Classification

Raster maps are geo-reference images, in which infomration is classified in the pixel values. In simple terms, the difference between regular images and geo-referenced raster images is that the latter has each pixel associated to a geographic (or projected) coordinate. So instead of just using rows and columns like we would do in regular images, we can access data in raster images using latitude and longitude.

The following example uses a land cover map for 2019 for Barton county, KS. This county is located in central Kansas and has a bit of everything: towns, lakes, roads, cropland, and grassland. The goal of this exercise is to learn how read raster data, access the individual values and coordinates, apply some operations to the data to extract features, and use image analysis to extract features of potential interest.

import numpy as np
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
import xarray as xr

# Load maps
B2010 = xr.open_rasterio('../datasets/land_cover/barton_county_ks_2010.tif')
B2019 = xr.open_rasterio('../datasets/land_cover/barton_county_ks_2019.tif')

# Show map
options = opts.Image(cmap='Spectral',aspect='equal',colorbar=True)
barton_2010_map = hv.Image( (B2010.x,B2010.y,B2010.values[0]),label='Barton County, KS 2010').opts(options)
barton_2019_map = hv.Image( (B2019.x,B2019.y,B2019.values[0]),label='Barton County, KS 2019' ).opts(options)

(barton_2019_map).opts(frame_width=300) # Use this to dsiplay both figure: (barton_2019_map + barton_2019_map).opts(frame_width=800).cols(2)

Map codes

The image has codes for each land cover encoded in uint8, which means unsigned integers of 8 bits. The values range from 0 to 255, meaning that this way of encoding information can only store 256 possible land covers, which is plenty for most applications. For example a pixel with a value of 111 representes areas occupied by water in 2019 and pixels with a value of 24 represents fields planted with winter wheat in 2019. The total of 256 land covers comes from the fact that a two-bit system would allow us to represent only two land covers using either 0 or 1. We actually use this format quite often when we plot boolean arrays.

To expand the number of possibilities we need more slots with 0s and 1s. So, since \(2^8=256\), using eight digits that can be either 0s or 1s allows for a total of 256 combinations. So, the value for grassland 176 would be represented in binary as 10110000 You can try this on your own by running the following Python code: bin(176). The first two values of the printed sequence (something like 0b) simply mean to tell you that the numbers are represented using the binary system.

Land Cover Codes

Land cover

Code

Comments

Background

0

Corn

1

Cotton

2

Rice

3

Sorghum

4

Soybeans

5

Wheat

24

Winter wheat and Spring wheat

Open water

111

Lakes, ponds, rivers, streams and creeks

Urban areas

121, 122, 123, 124

Each code is for very low to high development area

Barren land

131

Forests

141, 142, 143

Refers to deciduous, evergreen and wetland forests

Shrublands

152

Grassland/Pastures

176

Wetlands

190

 

# Print unique map codes
np.unique(B2019.values)

array([  0,   1,   2,   4,   5,   6,  21,  24,  26,  27,  28,  29,  31,
        36,  37,  44,  53,  58,  59,  61, 111, 121, 122, 123, 124, 131,
       141, 142, 143, 152, 176, 190, 195, 205, 225, 228, 236, 238],
      dtype=uint8)

# Classify some land covers
water_2019 = B2019 == 111
urban_2019 = (B2019>=121) & (B2019<=124)
wheat_2019 = B2019 == 24
grassland_2019 = B2019 == 176

# Plot land covers
water_map_2019 = hv.Image( (water_2019.x,water_2019.y,water_2019.values[0]) ).opts(cmap='Blues',title='Water bodies')
urban_map_2019 = hv.Image( (urban_2019.x,urban_2019.y,urban_2019.values[0]) ).opts(cmap='Purples',title='Urban area')
wheat_map_2019 = hv.Image( (wheat_2019.x,wheat_2019.y,wheat_2019.values[0]) ).opts(cmap='Reds',title='Wheat')
grassland_map_2019 = hv.Image( (grassland_2019.x,grassland_2019.y,grassland_2019.values[0]) ).opts(cmap='Greens',title='Grassland/Pasture')

(water_map_2019 + urban_map_2019  + wheat_map_2019 + grassland_map_2019).cols(2)

# Find area of each land cover
cell_area = 30*30 # m
water_area_2019 = water_2019.values.sum()*cell_area/10**6 # Area in km^2
water_area_2019 

36.7911

from skimage.measure import label, regionprops, regionprops_table
import pandas as pd

label_img = label(water_2019.values)
regions = regionprops(label_img)

props = regionprops_table(label_img[0], properties=('area','perimeter'))

props = pd.DataFrame(props)
props.head()
area perimeter
0 2 0.000000
1 17 16.899495
2 6 5.000000
3 4 4.000000
4 6 5.207107 
# Compare number of elements in props and labeled image.
print(props.shape)
print(np.unique(label_img[0]))

# Find total number of connected water bodies
print('Total number of water bodies in Barton county, KS in 2019:',props.shape[0])

(2612, 2)
[   0    1    2 ... 2610 2611 2612]
Total number of water bodies in Barton county, KS in 2019: 2612

# Plot all connected water bodies
hv.Image( (grassland_2019.x,grassland_2019.y,label_img[0]) ).opts(cmap='Colorblind_r', 
                                                                  frame_width=250,
                                                                  aspect='equal')

Practice

  • Quantify the decadal change in urban area and lake size using the land cover maps for 2010 vs 2019. Express your answer in \(km^2\).

  • Create a map that shows the change by the lake

  • Quantify the decadal change in the area of total cropland (wheat,corn, soybeans) and grassland. Express your answer in \(km^2\).

Rasters from NASA Worldview Snapshots Denoise Yield Monitor Data